Spatial Stats Part 1 new updated version
Kostas Goulias
2/13/2018
This is a shortened version of the lab instructions eliminating a few statwide analyses
Part 1 Preliminary Work Setting Up Analysis of Census Block Groups
A few important libraries for spatial analysis are:
install.packages(“maptools”, dependencies = TRUE) install.packages(“spdep”, dependencies = TRUE) install.packages(“leaflet”, dependencies = TRUE) install.packages(“RColorBrewer”, dependencies = TRUE)
library(spdep)## Loading required package: sp## Loading required package: Matrix## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge')`library(maptools)## Checking rgeos availability: TRUElibrary(leaflet)
library(RColorBrewer)Read the data by setting the directory to point at the right place You will need to download a shp file from Gauchospace and save it in this directory
CA.poly <- readShapePoly('./shp/LPA_Pop_Char_bg.shp')## Warning: use rgdal::readOGR or sf::st_readWe use the function class to verify the we have a spatial polygons data frame
class(CA.poly)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"A SpatialPolygonsDataFrame object brings together the spatial representations of the polygons with data. data: Object of class “data.frame”; attribute table polygons: Object of class “list”; see SpatialPolygons-class plotOrder: Object of class “integer”; see SpatialPolygons-class bbox: Object of class “matrix”; see Spatial-class proj4string: Object of class “CRS”; see CRS-class
The identifying tags of the polygons in the slot are matched with the row names of the data frame to make sure that the correct data rows are associated with the correct spatial object. In this example we have US Census block groups and for each block group we have added the behvior of people and the chractersitics of the households and persons
A Spatial Polygon Data Frame has four components (in the R jargon these called slots)
Component 1: Data contains the variables that are used in the analysis such as number of households with zero cars, vehicle miles of travel. remember these are at the block group level
str(slot(CA.poly, "data"))## 'data.frame': 23198 obs. of 105 variables:
## $ OBJECTID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ STATEFP : Factor w/ 1 level "06": 1 1 1 1 1 1 1 1 1 1 ...
## $ COUNTYFP : Factor w/ 58 levels "001","003","005",..: 42 37 30 12 19 30 30 19 19 33 ...
## $ TRACTCE : Factor w/ 6521 levels "000100","000101",..: 319 1482 2412 1083 5937 2492 1898 3170 3757 2302 ...
## $ BLKGRPCE : Factor w/ 10 levels "0","1","2","3",..: 2 3 3 2 2 3 6 2 2 2 ...
## $ AFFGEOID : Factor w/ 23198 levels "1500000US060014001001",..: 19731 17394 11471 2765 8865 11726 11313 4165 5247 13779 ...
## $ NAME : Factor w/ 10 levels "0","1","2","3",..: 2 3 3 2 2 3 6 2 2 2 ...
## $ LSAD : Factor w/ 1 level "BG": 1 1 1 1 1 1 1 1 1 1 ...
## $ ALAND : num 6.70e+06 3.80e+06 2.89e+05 8.03e+08 2.61e+05 ...
## $ AWATER : num 502651 2039168 0 1577551 2184908 ...
## $ GEO_ID2 : num 6.08e+10 6.07e+10 6.06e+10 6.02e+10 6.04e+10 ...
## $ GEOID_1 : num 6.08e+10 6.07e+10 6.06e+10 6.02e+10 6.04e+10 ...
## $ HHAGE1 : num 3 14 0 21 50 5 7 4 6 11 ...
## $ HHAGE2 : num 12 87 0 81 206 82 105 20 31 45 ...
## $ HHAGE3 : num 34 180 0 67 117 73 193 63 51 57 ...
## $ HHAGE4 : num 103 225 6 113 80 96 75 75 49 82 ...
## $ HHAGE5 : num 67 96 20 68 27 59 22 52 30 50 ...
## $ HHAGE6 : num 57 97 18 79 22 43 27 40 11 85 ...
## $ HHAGE7 : num 99 119 107 59 16 45 15 39 18 256 ...
## $ HHAGE8 : num 60 67 221 48 7 25 8 16 7 187 ...
## $ HHAGE9 : num 28 39 105 7 4 9 1 16 2 52 ...
## $ HHCHILD1 : num 336 643 477 431 477 409 199 210 92 706 ...
## $ HHCHILD2 : num 127 281 0 112 52 28 254 115 113 119 ...
## $ HHSIZE1 : num 80 219 250 196 322 238 80 68 18 415 ...
## $ HHSIZE2 : num 218 343 210 214 155 165 107 90 31 271 ...
## $ HHSIZE3 : num 67 153 17 64 31 25 100 55 32 68 ...
## $ HHSIZE4 : num 47 144 0 38 15 9 110 62 35 32 ...
## $ HHSIZE5 : num 33 50 0 24 6 0 44 26 30 23 ...
## $ HHSIZE6 : num 10 15 0 3 0 0 10 12 15 12 ...
## $ HHSIZE7 : num 8 0 0 4 0 0 2 12 44 4 ...
## $ rHHINC1 : num 11 60 0 82 11 17 0 0 0 37 ...
## $ rHHINC2 : num 11 129 62 136 76 58 0 0 50 255 ...
## $ rHHINC3 : num 50 39 130 144 114 23 0 49 32 47 ...
## $ rHHINC4 : num 16 67 179 105 168 6 12 17 39 178 ...
## $ rHHINC5 : num 36 69 42 76 48 101 12 95 26 114 ...
## $ rHHINC6 : num 61 173 47 0 63 50 73 67 27 55 ...
## $ rHHINC7 : num 79 121 9 0 25 40 173 60 31 34 ...
## $ rHHINC8 : num 67 75 0 0 0 81 67 0 0 39 ...
## $ rHHINC9 : num 132 191 8 0 24 61 116 37 0 66 ...
## $ n_hh : num 463 924 477 543 529 437 453 325 205 825 ...
## $ AGE1 : num 163 408 0 138 63 33 447 172 227 153 ...
## $ AGE2 : num 139 174 1 93 86 25 73 109 172 101 ...
## $ AGE3 : num 88 409 1 202 390 211 436 155 206 153 ...
## $ AGE4 : num 136 385 9 140 119 110 216 158 127 116 ...
## $ AGE5 : num 350 524 99 381 116 190 117 224 101 310 ...
## $ AGE6 : num 315 380 611 180 41 110 40 135 52 673 ...
## $ SEX1 : num 563 1077 258 495 402 ...
## $ SEX2 : num 628 1203 463 639 413 ...
## $ n_pr : num 1191 2280 721 1134 815 ...
## $ Shape_Leng: num 0.1822 0.1509 0.023 2.9636 0.0367 ...
## $ Shape_Area: num 6.64e-04 3.92e-04 2.81e-05 8.61e-02 2.33e-05 ...
## $ FID_ : int 0 1 2 3 4 5 6 7 8 9 ...
## $ STATEFP_1 : Factor w/ 1 level "06": 1 1 1 1 1 1 1 1 1 1 ...
## $ COUNTYFP_1: Factor w/ 58 levels "001","003","005",..: 42 37 30 12 19 30 30 19 19 33 ...
## $ TRACTCE_1 : Factor w/ 6521 levels "000100","000101",..: 319 1482 2412 1083 5937 2492 1898 3170 3757 2302 ...
## $ BLKGRPCE_1: Factor w/ 10 levels "0","1","2","3",..: 2 3 3 2 2 3 6 2 2 2 ...
## $ AFFGEOID_1: Factor w/ 23198 levels "1500000US060014001001",..: 19731 17394 11471 2765 8865 11726 11313 4165 5247 13779 ...
## $ ALAND_1 : num 6.70e+06 3.80e+06 2.89e+05 8.03e+08 2.61e+05 ...
## $ AWATER_1 : num 502651 2039168 0 1577551 2184908 ...
## $ Shape_Ar_1: num 6.64e-04 3.92e-04 2.81e-05 8.61e-02 2.33e-05 ...
## $ GEO_ID2_1 : num 6.08e+10 6.07e+10 6.06e+10 6.02e+10 6.04e+10 ...
## $ pumano : int 8303 7304 5904 2300 3765 5903 5915 3724 3745 6515 ...
## $ countyname: Factor w/ 58 levels "Alameda","Alpine",..: 42 37 30 12 19 30 30 19 19 33 ...
## $ GEOID : num 6.08e+10 6.07e+10 6.06e+10 6.02e+10 6.04e+10 ...
## $ n_hh_1 : num 463 924 477 543 529 437 453 325 205 825 ...
## $ n_pr_1 : num 1191 2280 721 1134 815 ...
## $ walkM : num 181 452 67 160 418 ...
## $ bikeM : num 232.6 186.2 34 37.5 384.9 ...
## $ drvalM : num 14067 25300 5429 11066 9466 ...
## $ drvothM : num 6583 13495 4246 8737 3667 ...
## $ passM : num 7720 19778 3538 10707 4177 ...
## $ planeM : num 922 838 0 335 0 ...
## $ otherM : num 643 4394 439 598 2144 ...
## $ allM : num 30348 64443 13754 31640 20256 ...
## $ walkD : num 3972 11498 2992 2865 10470 ...
## $ bikeD : num 1812 2448 304 289 2910 ...
## $ drvalD : num 34663 60042 14727 23800 28254 ...
## $ drvothD : num 14518 28376 8600 17921 8804 ...
## $ passD : num 17764 40452 7514 22878 10399 ...
## $ planeD : num 1577 4256 0 1200 0 ...
## $ otherD : num 1917 9406 1195 1946 6084 ...
## $ allD : num 76223 156478 35332 70899 66921 ...
## $ walkT : num 281 879 174 180 848 388 498 287 566 566 ...
## $ bikeT : num 44 135 10 10 126 18 46 48 17 59 ...
## $ drvalT : num 1725 2887 905 1099 1379 ...
## $ drvothT : num 687 1347 342 680 407 ...
## $ passT : num 865 1780 299 882 459 254 1990 761 833 889 ...
## $ planeT : num 11 11 0 4 0 3 18 4 1 11 ...
## $ otherT : num 67 282 52 50 249 171 136 72 99 166 ...
## $ allT : num 3680 7321 1782 2905 3468 ...
## $ VMT_drv : num 20650 38795 9675 19802 13132 ...
## $ VMT_pas : num 7720 19778 3538 10707 4177 ...
## $ VMT : num 28370 58573 13214 30509 17309 ...
## $ pr_LDTrips: num 86 170 60 133 39 74 197 131 26 155 ...
## $ HHVEH0 : num 6 61 33 24 91 34 5 4 18 92 ...
## $ HHVEH1 : num 102 264 300 248 310 210 94 88 81 413 ...
## $ HHVEH2 : num 244 411 130 194 114 159 233 152 53 246 ...
## $ HHVEH3 : num 87 136 11 57 14 23 93 62 34 50 ...
## $ HHVEH4 : num 17 41 1 16 0 7 23 16 13 17 ...
## [list output truncated]
## - attr(*, "data_types")= chr "N" "C" "C" "C" ...summary(slot(CA.poly, "data"))## OBJECTID STATEFP COUNTYFP TRACTCE BLKGRPCE
## Min. : 1 06:23198 037 :6423 000200 : 71 1 :8036
## 1st Qu.: 5800 059 :1822 000300 : 70 2 :7353
## Median :11600 073 :1794 000400 : 68 3 :4736
## Mean :11600 071 :1092 000600 : 56 4 :2139
## 3rd Qu.:17399 085 :1075 000800 : 55 5 : 710
## Max. :23198 001 :1047 000100 : 51 6 : 169
## (Other):9945 (Other):22827 (Other): 55
## AFFGEOID NAME LSAD
## 1500000US060014001001: 1 1 :8036 BG:23198
## 1500000US060014002001: 1 2 :7353
## 1500000US060014002002: 1 3 :4736
## 1500000US060014003001: 1 4 :2139
## 1500000US060014003002: 1 5 : 710
## 1500000US060014003003: 1 6 : 169
## (Other) :23192 (Other): 55
## ALAND AWATER GEO_ID2
## Min. :0.000e+00 Min. :0.000e+00 Min. :6.001e+10
## 1st Qu.:2.988e+05 1st Qu.:0.000e+00 1st Qu.:6.037e+10
## Median :5.461e+05 Median :0.000e+00 Median :6.059e+10
## Mean :1.739e+07 Mean :5.722e+05 Mean :6.055e+10
## 3rd Qu.:1.296e+06 3rd Qu.:0.000e+00 3rd Qu.:6.073e+10
## Max. :1.610e+10 Max. :2.406e+09 Max. :6.115e+10
##
## GEOID_1 HHAGE1 HHAGE2 HHAGE3
## Min. :0.000e+00 Min. : 0.00 Min. : 0.00 Min. : 0.0
## 1st Qu.:6.037e+10 1st Qu.: 5.00 1st Qu.: 36.00 1st Qu.: 63.0
## Median :6.059e+10 Median : 12.00 Median : 63.00 Median : 91.0
## Mean :6.039e+10 Mean : 21.89 Mean : 85.11 Mean : 107.9
## 3rd Qu.:6.073e+10 3rd Qu.: 25.00 3rd Qu.: 107.00 3rd Qu.: 133.0
## Max. :6.115e+10 Max. :2130.00 Max. :2879.00 Max. :1878.0
##
## HHAGE4 HHAGE5 HHAGE6 HHAGE7
## Min. : 0.0 Min. : 0.00 Min. : 0.0 Min. : 0.0
## 1st Qu.: 76.0 1st Qu.: 33.00 1st Qu.: 27.0 1st Qu.: 31.0
## Median : 105.0 Median : 46.00 Median : 39.0 Median : 48.0
## Mean : 119.6 Mean : 52.21 Mean : 44.8 Mean : 57.3
## 3rd Qu.: 147.0 3rd Qu.: 65.00 3rd Qu.: 56.0 3rd Qu.: 72.0
## Max. :1361.0 Max. :381.00 Max. :454.0 Max. :828.0
##
## HHAGE8 HHAGE9 HHCHILD1 HHCHILD2
## Min. : 0.0 Min. : 0.00 Min. : 0.0 Min. : 0
## 1st Qu.: 17.0 1st Qu.: 5.00 1st Qu.: 197.0 1st Qu.: 116
## Median : 28.0 Median : 10.00 Median : 291.0 Median : 174
## Mean : 36.8 Mean : 16.48 Mean : 339.1 Mean : 203
## 3rd Qu.: 46.0 3rd Qu.: 19.00 3rd Qu.: 428.0 3rd Qu.: 254
## Max. :781.0 Max. :510.00 Max. :3363.0 Max. :4664
##
## HHSIZE1 HHSIZE2 HHSIZE3 HHSIZE4
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.00
## 1st Qu.: 54.0 1st Qu.: 88.0 1st Qu.: 54.0 1st Qu.: 47.00
## Median : 92.0 Median : 136.0 Median : 77.0 Median : 69.00
## Mean : 126.1 Mean : 157.5 Mean : 88.2 Mean : 81.28
## 3rd Qu.: 159.0 3rd Qu.: 202.0 3rd Qu.: 110.0 3rd Qu.: 101.00
## Max. :1884.0 Max. :1519.0 Max. :1678.0 Max. :1770.00
##
## HHSIZE5 HHSIZE6 HHSIZE7 rHHINC1
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 23.00 1st Qu.: 8.00 1st Qu.: 4.00 1st Qu.: 0.00
## Median : 38.00 Median : 17.00 Median : 12.00 Median : 17.00
## Mean : 44.93 Mean : 22.01 Mean : 22.11 Mean : 28.75
## 3rd Qu.: 58.00 3rd Qu.: 30.00 3rd Qu.: 32.00 3rd Qu.: 41.00
## Max. :830.00 Max. :412.00 Max. :418.00 Max. :722.00
##
## rHHINC2 rHHINC3 rHHINC4 rHHINC5
## Min. : 0.00 Min. : 0.0 Min. : 0.00 Min. : 0.00
## 1st Qu.: 26.00 1st Qu.: 16.0 1st Qu.: 29.00 1st Qu.: 47.00
## Median : 59.00 Median : 39.0 Median : 56.00 Median : 81.00
## Mean : 78.56 Mean : 49.3 Mean : 68.37 Mean : 95.38
## 3rd Qu.: 110.00 3rd Qu.: 71.0 3rd Qu.: 95.00 3rd Qu.: 127.00
## Max. :1018.00 Max. :713.0 Max. :1598.00 Max. :1611.00
##
## rHHINC6 rHHINC7 rHHINC8 rHHINC9
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 29.00 1st Qu.: 28.00 1st Qu.: 0.00 1st Qu.: 0.00
## Median : 56.00 Median : 64.00 Median : 19.00 Median : 11.00
## Mean : 69.41 Mean : 81.93 Mean : 34.87 Mean : 35.54
## 3rd Qu.: 95.00 3rd Qu.: 114.00 3rd Qu.: 50.00 3rd Qu.: 45.00
## Max. :1201.00 Max. :1415.00 Max. :668.00 Max. :1270.00
##
## n_hh AGE1 AGE2 AGE3
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 350.0 1st Qu.: 166.0 1st Qu.: 122.0 1st Qu.: 181.0
## Median : 480.0 Median : 261.0 Median : 193.0 Median : 276.0
## Mean : 542.1 Mean : 315.4 Mean : 224.7 Mean : 326.9
## 3rd Qu.: 665.0 3rd Qu.: 403.0 3rd Qu.: 288.0 3rd Qu.: 410.0
## Max. :6016.0 Max. :4934.0 Max. :9156.0 Max. :6244.0
##
## AGE4 AGE5 AGE6 SEX1
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 138.0 1st Qu.: 180.0 1st Qu.: 100.0 1st Qu.: 480.0
## Median : 194.0 Median : 247.0 Median : 151.0 Median : 666.0
## Mean : 221.6 Mean : 279.7 Mean : 184.1 Mean : 748.9
## 3rd Qu.: 271.0 3rd Qu.: 343.0 3rd Qu.: 230.0 3rd Qu.: 919.0
## Max. :3270.0 Max. :2320.0 Max. :2689.0 Max. :11511.0
##
## SEX2 n_pr Shape_Leng Shape_Area
## Min. : 0.0 Min. : 0 Min. :0.000039 Min. :0.0000000
## 1st Qu.: 518.0 1st Qu.: 1001 1st Qu.:0.024613 1st Qu.:0.0000297
## Median : 715.0 Median : 1380 Median :0.034106 Median :0.0000549
## Mean : 803.5 Mean : 1552 Mean :0.088440 Mean :0.0017962
## 3rd Qu.: 980.0 3rd Qu.: 1895 3rd Qu.:0.058370 3rd Qu.:0.0001306
## Max. :10054.0 Max. :20885 Max. :8.865723 Max. :1.5939843
##
## FID_ STATEFP_1 COUNTYFP_1 TRACTCE_1 BLKGRPCE_1
## Min. : 0 06:23198 037 :6423 000200 : 71 1 :8036
## 1st Qu.: 5799 059 :1822 000300 : 70 2 :7353
## Median :11598 073 :1794 000400 : 68 3 :4736
## Mean :11598 071 :1092 000600 : 56 4 :2139
## 3rd Qu.:17398 085 :1075 000800 : 55 5 : 710
## Max. :23197 001 :1047 000100 : 51 6 : 169
## (Other):9945 (Other):22827 (Other): 55
## AFFGEOID_1 ALAND_1 AWATER_1
## 1500000US060014001001: 1 Min. :0.000e+00 Min. :0.000e+00
## 1500000US060014002001: 1 1st Qu.:2.988e+05 1st Qu.:0.000e+00
## 1500000US060014002002: 1 Median :5.461e+05 Median :0.000e+00
## 1500000US060014003001: 1 Mean :1.739e+07 Mean :5.722e+05
## 1500000US060014003002: 1 3rd Qu.:1.296e+06 3rd Qu.:0.000e+00
## 1500000US060014003003: 1 Max. :1.610e+10 Max. :2.406e+09
## (Other) :23192
## Shape_Ar_1 GEO_ID2_1 pumano
## Min. :0.0000000 Min. :6.001e+10 Min. : 0
## 1st Qu.:0.0000297 1st Qu.:6.037e+10 1st Qu.: 3723
## Median :0.0000549 Median :6.059e+10 Median : 5905
## Mean :0.0017962 Mean :6.055e+10 Mean : 5425
## 3rd Qu.:0.0001306 3rd Qu.:6.073e+10 3rd Qu.: 7316
## Max. :1.5939843 Max. :6.115e+10 Max. :11300
##
## countyname GEOID n_hh_1 n_pr_1
## Los Angeles :6423 Min. :0.000e+00 Min. : 0.0 Min. : 0
## Orange :1822 1st Qu.:6.037e+10 1st Qu.: 350.0 1st Qu.: 1001
## San Diego :1794 Median :6.059e+10 Median : 480.0 Median : 1380
## San Bernadino:1092 Mean :6.039e+10 Mean : 542.1 Mean : 1552
## Santa Clara :1075 3rd Qu.:6.073e+10 3rd Qu.: 665.0 3rd Qu.: 1895
## Alameda :1047 Max. :6.115e+10 Max. :6016.0 Max. :20885
## (Other) :9945
## walkM bikeM drvalM drvothM
## Min. : 0.0 Min. : 0.00 Min. : 0 Min. : 0
## 1st Qu.: 206.7 1st Qu.: 72.41 1st Qu.: 8560 1st Qu.: 5180
## Median : 330.6 Median : 127.06 Median : 12131 Median : 7482
## Mean : 398.1 Mean : 163.46 Mean : 14282 Mean : 8926
## 3rd Qu.: 511.1 3rd Qu.: 210.11 3rd Qu.: 17503 3rd Qu.: 10852
## Max. :3697.9 Max. :4100.94 Max. :205090 Max. :266464
##
## passM planeM otherM allM
## Min. : 0 Min. : 0.00 Min. : 0 Min. : 0
## 1st Qu.: 7031 1st Qu.: 83.78 1st Qu.: 1157 1st Qu.: 24199
## Median : 10824 Median : 335.11 Median : 1819 Median : 34459
## Mean : 14608 Mean : 633.93 Mean : 2168 Mean : 41180
## 3rd Qu.: 17316 3rd Qu.: 837.77 3rd Qu.: 2744 3rd Qu.: 49962
## Max. :676613 Max. :16839.22 Max. :38781 Max. :1167862
##
## walkD bikeD drvalD drvothD
## Min. : 0 Min. : 0 Min. : 0 Min. : 0
## 1st Qu.: 5017 1st Qu.: 689 1st Qu.: 21739 1st Qu.: 11866
## Median : 8116 Median : 1169 Median : 30139 Median : 16733
## Mean : 9850 Mean : 1496 Mean : 34575 Mean : 19379
## 3rd Qu.: 12692 3rd Qu.: 1912 3rd Qu.: 42305 3rd Qu.: 23536
## Max. :128210 Max. :29851 Max. :421362 Max. :314652
##
## passD planeD otherD allD
## Min. : 0 Min. : 0 Min. : 0 Min. : 0
## 1st Qu.: 17602 1st Qu.: 170 1st Qu.: 3692 1st Qu.: 66732
## Median : 26111 Median : 805 Median : 5760 Median : 92364
## Mean : 31697 Mean : 1426 Mean : 6797 Mean : 105219
## 3rd Qu.: 38562 3rd Qu.: 2034 3rd Qu.: 8756 3rd Qu.: 127822
## Max. :1035072 Max. :25910 Max. :72248 Max. :1897854
##
## walkT bikeT drvalT drvothT
## Min. : 0.0 Min. : 0.00 Min. : 0 Min. : 0.0
## 1st Qu.: 391.0 1st Qu.: 34.00 1st Qu.: 1036 1st Qu.: 573.0
## Median : 646.5 Median : 60.00 Median : 1435 Median : 802.0
## Mean : 798.4 Mean : 76.72 Mean : 1639 Mean : 919.3
## 3rd Qu.:1038.0 3rd Qu.: 100.00 3rd Qu.: 2011 3rd Qu.: 1126.0
## Max. :8072.0 Max. :1461.00 Max. :19585 Max. :11940.0
##
## passT planeT otherT allT
## Min. : 0 Min. : 0.00 Min. : 0.0 Min. : 0
## 1st Qu.: 799 1st Qu.: 1.00 1st Qu.: 128.0 1st Qu.: 3240
## Median : 1159 Median : 4.00 Median : 208.0 Median : 4482
## Mean : 1352 Mean : 7.42 Mean : 256.1 Mean : 5048
## 3rd Qu.: 1665 3rd Qu.: 10.00 3rd Qu.: 335.0 3rd Qu.: 6192
## Max. :18865 Max. :201.00 Max. :2975.0 Max. :59499
##
## VMT_drv VMT_pas VMT pr_LDTrips
## Min. : 0 Min. : 0 Min. : 0 Min. : 0.0
## 1st Qu.: 14022 1st Qu.: 7031 1st Qu.: 21738 1st Qu.: 67.0
## Median : 19697 Median : 10824 Median : 31289 Median : 105.0
## Mean : 23208 Mean : 14608 Mean : 37816 Mean : 131.6
## 3rd Qu.: 28201 3rd Qu.: 17316 3rd Qu.: 45815 3rd Qu.: 164.0
## Max. :471555 Max. :676613 Max. :1148168 Max. :5039.0
##
## HHVEH0 HHVEH1 HHVEH2 HHVEH3
## Min. : 0.00 Min. : 0.0 Min. : 0.0 Min. : 0.00
## 1st Qu.: 12.00 1st Qu.: 102.0 1st Qu.: 131.0 1st Qu.: 42.00
## Median : 27.00 Median : 153.0 Median : 183.0 Median : 63.00
## Mean : 39.65 Mean : 183.8 Mean : 212.3 Mean : 76.62
## 3rd Qu.: 53.00 3rd Qu.: 233.0 3rd Qu.: 260.0 3rd Qu.: 94.00
## Max. :760.00 Max. :1824.0 Max. :2997.0 Max. :1620.00
##
## HHVEH4 HHVEH5 HHVEH6 HHVEH7
## Min. : 0.00 Min. : 0.00 Min. : 0.000 Min. : 0.0000
## 1st Qu.: 10.00 1st Qu.: 1.00 1st Qu.: 0.000 1st Qu.: 0.0000
## Median : 17.00 Median : 3.00 Median : 1.000 Median : 0.0000
## Mean : 21.09 Mean : 4.95 Mean : 2.787 Mean : 0.7652
## 3rd Qu.: 27.00 3rd Qu.: 6.00 3rd Qu.: 3.000 3rd Qu.: 1.0000
## Max. :320.00 Max. :149.00 Max. :72.000 Max. :24.0000
##
## HHVEH8 notrav LPAgrp
## Min. : 0.0000 Min. : 0.0 Min. :0.000
## 1st Qu.: 0.0000 1st Qu.: 226.0 1st Qu.:3.000
## Median : 0.0000 Median : 321.0 Median :3.000
## Mean : 0.1916 Mean : 359.6 Mean :3.178
## 3rd Qu.: 0.0000 3rd Qu.: 447.0 3rd Qu.:4.000
## Max. :38.0000 Max. :3763.0 Max. :4.000
## This shows an output like : ‘data.frame’: 23198 obs. of 105 variables:
Component 2: This is the polygon slot and contains the “shape” information.
The following will create a map of all the polygons representing block groups
plot(CA.poly)
Analysis of the data in these blockgroups
summary(CA.poly)## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -124.40959 -114.13121
## y 32.53416 42.00952
## Is projected: NA
## proj4string : [NA]
## Data attributes:
## OBJECTID STATEFP COUNTYFP TRACTCE BLKGRPCE
## Min. : 1 06:23198 037 :6423 000200 : 71 1 :8036
## 1st Qu.: 5800 059 :1822 000300 : 70 2 :7353
## Median :11600 073 :1794 000400 : 68 3 :4736
## Mean :11600 071 :1092 000600 : 56 4 :2139
## 3rd Qu.:17399 085 :1075 000800 : 55 5 : 710
## Max. :23198 001 :1047 000100 : 51 6 : 169
## (Other):9945 (Other):22827 (Other): 55
## AFFGEOID NAME LSAD
## 1500000US060014001001: 1 1 :8036 BG:23198
## 1500000US060014002001: 1 2 :7353
## 1500000US060014002002: 1 3 :4736
## 1500000US060014003001: 1 4 :2139
## 1500000US060014003002: 1 5 : 710
## 1500000US060014003003: 1 6 : 169
## (Other) :23192 (Other): 55
## ALAND AWATER GEO_ID2
## Min. :0.000e+00 Min. :0.000e+00 Min. :6.001e+10
## 1st Qu.:2.988e+05 1st Qu.:0.000e+00 1st Qu.:6.037e+10
## Median :5.461e+05 Median :0.000e+00 Median :6.059e+10
## Mean :1.739e+07 Mean :5.722e+05 Mean :6.055e+10
## 3rd Qu.:1.296e+06 3rd Qu.:0.000e+00 3rd Qu.:6.073e+10
## Max. :1.610e+10 Max. :2.406e+09 Max. :6.115e+10
##
## GEOID_1 HHAGE1 HHAGE2 HHAGE3
## Min. :0.000e+00 Min. : 0.00 Min. : 0.00 Min. : 0.0
## 1st Qu.:6.037e+10 1st Qu.: 5.00 1st Qu.: 36.00 1st Qu.: 63.0
## Median :6.059e+10 Median : 12.00 Median : 63.00 Median : 91.0
## Mean :6.039e+10 Mean : 21.89 Mean : 85.11 Mean : 107.9
## 3rd Qu.:6.073e+10 3rd Qu.: 25.00 3rd Qu.: 107.00 3rd Qu.: 133.0
## Max. :6.115e+10 Max. :2130.00 Max. :2879.00 Max. :1878.0
##
## HHAGE4 HHAGE5 HHAGE6 HHAGE7
## Min. : 0.0 Min. : 0.00 Min. : 0.0 Min. : 0.0
## 1st Qu.: 76.0 1st Qu.: 33.00 1st Qu.: 27.0 1st Qu.: 31.0
## Median : 105.0 Median : 46.00 Median : 39.0 Median : 48.0
## Mean : 119.6 Mean : 52.21 Mean : 44.8 Mean : 57.3
## 3rd Qu.: 147.0 3rd Qu.: 65.00 3rd Qu.: 56.0 3rd Qu.: 72.0
## Max. :1361.0 Max. :381.00 Max. :454.0 Max. :828.0
##
## HHAGE8 HHAGE9 HHCHILD1 HHCHILD2
## Min. : 0.0 Min. : 0.00 Min. : 0.0 Min. : 0
## 1st Qu.: 17.0 1st Qu.: 5.00 1st Qu.: 197.0 1st Qu.: 116
## Median : 28.0 Median : 10.00 Median : 291.0 Median : 174
## Mean : 36.8 Mean : 16.48 Mean : 339.1 Mean : 203
## 3rd Qu.: 46.0 3rd Qu.: 19.00 3rd Qu.: 428.0 3rd Qu.: 254
## Max. :781.0 Max. :510.00 Max. :3363.0 Max. :4664
##
## HHSIZE1 HHSIZE2 HHSIZE3 HHSIZE4
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.00
## 1st Qu.: 54.0 1st Qu.: 88.0 1st Qu.: 54.0 1st Qu.: 47.00
## Median : 92.0 Median : 136.0 Median : 77.0 Median : 69.00
## Mean : 126.1 Mean : 157.5 Mean : 88.2 Mean : 81.28
## 3rd Qu.: 159.0 3rd Qu.: 202.0 3rd Qu.: 110.0 3rd Qu.: 101.00
## Max. :1884.0 Max. :1519.0 Max. :1678.0 Max. :1770.00
##
## HHSIZE5 HHSIZE6 HHSIZE7 rHHINC1
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 23.00 1st Qu.: 8.00 1st Qu.: 4.00 1st Qu.: 0.00
## Median : 38.00 Median : 17.00 Median : 12.00 Median : 17.00
## Mean : 44.93 Mean : 22.01 Mean : 22.11 Mean : 28.75
## 3rd Qu.: 58.00 3rd Qu.: 30.00 3rd Qu.: 32.00 3rd Qu.: 41.00
## Max. :830.00 Max. :412.00 Max. :418.00 Max. :722.00
##
## rHHINC2 rHHINC3 rHHINC4 rHHINC5
## Min. : 0.00 Min. : 0.0 Min. : 0.00 Min. : 0.00
## 1st Qu.: 26.00 1st Qu.: 16.0 1st Qu.: 29.00 1st Qu.: 47.00
## Median : 59.00 Median : 39.0 Median : 56.00 Median : 81.00
## Mean : 78.56 Mean : 49.3 Mean : 68.37 Mean : 95.38
## 3rd Qu.: 110.00 3rd Qu.: 71.0 3rd Qu.: 95.00 3rd Qu.: 127.00
## Max. :1018.00 Max. :713.0 Max. :1598.00 Max. :1611.00
##
## rHHINC6 rHHINC7 rHHINC8 rHHINC9
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 29.00 1st Qu.: 28.00 1st Qu.: 0.00 1st Qu.: 0.00
## Median : 56.00 Median : 64.00 Median : 19.00 Median : 11.00
## Mean : 69.41 Mean : 81.93 Mean : 34.87 Mean : 35.54
## 3rd Qu.: 95.00 3rd Qu.: 114.00 3rd Qu.: 50.00 3rd Qu.: 45.00
## Max. :1201.00 Max. :1415.00 Max. :668.00 Max. :1270.00
##
## n_hh AGE1 AGE2 AGE3
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 350.0 1st Qu.: 166.0 1st Qu.: 122.0 1st Qu.: 181.0
## Median : 480.0 Median : 261.0 Median : 193.0 Median : 276.0
## Mean : 542.1 Mean : 315.4 Mean : 224.7 Mean : 326.9
## 3rd Qu.: 665.0 3rd Qu.: 403.0 3rd Qu.: 288.0 3rd Qu.: 410.0
## Max. :6016.0 Max. :4934.0 Max. :9156.0 Max. :6244.0
##
## AGE4 AGE5 AGE6 SEX1
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 0.0
## 1st Qu.: 138.0 1st Qu.: 180.0 1st Qu.: 100.0 1st Qu.: 480.0
## Median : 194.0 Median : 247.0 Median : 151.0 Median : 666.0
## Mean : 221.6 Mean : 279.7 Mean : 184.1 Mean : 748.9
## 3rd Qu.: 271.0 3rd Qu.: 343.0 3rd Qu.: 230.0 3rd Qu.: 919.0
## Max. :3270.0 Max. :2320.0 Max. :2689.0 Max. :11511.0
##
## SEX2 n_pr Shape_Leng Shape_Area
## Min. : 0.0 Min. : 0 Min. :0.000039 Min. :0.0000000
## 1st Qu.: 518.0 1st Qu.: 1001 1st Qu.:0.024613 1st Qu.:0.0000297
## Median : 715.0 Median : 1380 Median :0.034106 Median :0.0000549
## Mean : 803.5 Mean : 1552 Mean :0.088440 Mean :0.0017962
## 3rd Qu.: 980.0 3rd Qu.: 1895 3rd Qu.:0.058370 3rd Qu.:0.0001306
## Max. :10054.0 Max. :20885 Max. :8.865723 Max. :1.5939843
##
## FID_ STATEFP_1 COUNTYFP_1 TRACTCE_1 BLKGRPCE_1
## Min. : 0 06:23198 037 :6423 000200 : 71 1 :8036
## 1st Qu.: 5799 059 :1822 000300 : 70 2 :7353
## Median :11598 073 :1794 000400 : 68 3 :4736
## Mean :11598 071 :1092 000600 : 56 4 :2139
## 3rd Qu.:17398 085 :1075 000800 : 55 5 : 710
## Max. :23197 001 :1047 000100 : 51 6 : 169
## (Other):9945 (Other):22827 (Other): 55
## AFFGEOID_1 ALAND_1 AWATER_1
## 1500000US060014001001: 1 Min. :0.000e+00 Min. :0.000e+00
## 1500000US060014002001: 1 1st Qu.:2.988e+05 1st Qu.:0.000e+00
## 1500000US060014002002: 1 Median :5.461e+05 Median :0.000e+00
## 1500000US060014003001: 1 Mean :1.739e+07 Mean :5.722e+05
## 1500000US060014003002: 1 3rd Qu.:1.296e+06 3rd Qu.:0.000e+00
## 1500000US060014003003: 1 Max. :1.610e+10 Max. :2.406e+09
## (Other) :23192
## Shape_Ar_1 GEO_ID2_1 pumano
## Min. :0.0000000 Min. :6.001e+10 Min. : 0
## 1st Qu.:0.0000297 1st Qu.:6.037e+10 1st Qu.: 3723
## Median :0.0000549 Median :6.059e+10 Median : 5905
## Mean :0.0017962 Mean :6.055e+10 Mean : 5425
## 3rd Qu.:0.0001306 3rd Qu.:6.073e+10 3rd Qu.: 7316
## Max. :1.5939843 Max. :6.115e+10 Max. :11300
##
## countyname GEOID n_hh_1 n_pr_1
## Los Angeles :6423 Min. :0.000e+00 Min. : 0.0 Min. : 0
## Orange :1822 1st Qu.:6.037e+10 1st Qu.: 350.0 1st Qu.: 1001
## San Diego :1794 Median :6.059e+10 Median : 480.0 Median : 1380
## San Bernadino:1092 Mean :6.039e+10 Mean : 542.1 Mean : 1552
## Santa Clara :1075 3rd Qu.:6.073e+10 3rd Qu.: 665.0 3rd Qu.: 1895
## Alameda :1047 Max. :6.115e+10 Max. :6016.0 Max. :20885
## (Other) :9945
## walkM bikeM drvalM drvothM
## Min. : 0.0 Min. : 0.00 Min. : 0 Min. : 0
## 1st Qu.: 206.7 1st Qu.: 72.41 1st Qu.: 8560 1st Qu.: 5180
## Median : 330.6 Median : 127.06 Median : 12131 Median : 7482
## Mean : 398.1 Mean : 163.46 Mean : 14282 Mean : 8926
## 3rd Qu.: 511.1 3rd Qu.: 210.11 3rd Qu.: 17503 3rd Qu.: 10852
## Max. :3697.9 Max. :4100.94 Max. :205090 Max. :266464
##
## passM planeM otherM allM
## Min. : 0 Min. : 0.00 Min. : 0 Min. : 0
## 1st Qu.: 7031 1st Qu.: 83.78 1st Qu.: 1157 1st Qu.: 24199
## Median : 10824 Median : 335.11 Median : 1819 Median : 34459
## Mean : 14608 Mean : 633.93 Mean : 2168 Mean : 41180
## 3rd Qu.: 17316 3rd Qu.: 837.77 3rd Qu.: 2744 3rd Qu.: 49962
## Max. :676613 Max. :16839.22 Max. :38781 Max. :1167862
##
## walkD bikeD drvalD drvothD
## Min. : 0 Min. : 0 Min. : 0 Min. : 0
## 1st Qu.: 5017 1st Qu.: 689 1st Qu.: 21739 1st Qu.: 11866
## Median : 8116 Median : 1169 Median : 30139 Median : 16733
## Mean : 9850 Mean : 1496 Mean : 34575 Mean : 19379
## 3rd Qu.: 12692 3rd Qu.: 1912 3rd Qu.: 42305 3rd Qu.: 23536
## Max. :128210 Max. :29851 Max. :421362 Max. :314652
##
## passD planeD otherD allD
## Min. : 0 Min. : 0 Min. : 0 Min. : 0
## 1st Qu.: 17602 1st Qu.: 170 1st Qu.: 3692 1st Qu.: 66732
## Median : 26111 Median : 805 Median : 5760 Median : 92364
## Mean : 31697 Mean : 1426 Mean : 6797 Mean : 105219
## 3rd Qu.: 38562 3rd Qu.: 2034 3rd Qu.: 8756 3rd Qu.: 127822
## Max. :1035072 Max. :25910 Max. :72248 Max. :1897854
##
## walkT bikeT drvalT drvothT
## Min. : 0.0 Min. : 0.00 Min. : 0 Min. : 0.0
## 1st Qu.: 391.0 1st Qu.: 34.00 1st Qu.: 1036 1st Qu.: 573.0
## Median : 646.5 Median : 60.00 Median : 1435 Median : 802.0
## Mean : 798.4 Mean : 76.72 Mean : 1639 Mean : 919.3
## 3rd Qu.:1038.0 3rd Qu.: 100.00 3rd Qu.: 2011 3rd Qu.: 1126.0
## Max. :8072.0 Max. :1461.00 Max. :19585 Max. :11940.0
##
## passT planeT otherT allT
## Min. : 0 Min. : 0.00 Min. : 0.0 Min. : 0
## 1st Qu.: 799 1st Qu.: 1.00 1st Qu.: 128.0 1st Qu.: 3240
## Median : 1159 Median : 4.00 Median : 208.0 Median : 4482
## Mean : 1352 Mean : 7.42 Mean : 256.1 Mean : 5048
## 3rd Qu.: 1665 3rd Qu.: 10.00 3rd Qu.: 335.0 3rd Qu.: 6192
## Max. :18865 Max. :201.00 Max. :2975.0 Max. :59499
##
## VMT_drv VMT_pas VMT pr_LDTrips
## Min. : 0 Min. : 0 Min. : 0 Min. : 0.0
## 1st Qu.: 14022 1st Qu.: 7031 1st Qu.: 21738 1st Qu.: 67.0
## Median : 19697 Median : 10824 Median : 31289 Median : 105.0
## Mean : 23208 Mean : 14608 Mean : 37816 Mean : 131.6
## 3rd Qu.: 28201 3rd Qu.: 17316 3rd Qu.: 45815 3rd Qu.: 164.0
## Max. :471555 Max. :676613 Max. :1148168 Max. :5039.0
##
## HHVEH0 HHVEH1 HHVEH2 HHVEH3
## Min. : 0.00 Min. : 0.0 Min. : 0.0 Min. : 0.00
## 1st Qu.: 12.00 1st Qu.: 102.0 1st Qu.: 131.0 1st Qu.: 42.00
## Median : 27.00 Median : 153.0 Median : 183.0 Median : 63.00
## Mean : 39.65 Mean : 183.8 Mean : 212.3 Mean : 76.62
## 3rd Qu.: 53.00 3rd Qu.: 233.0 3rd Qu.: 260.0 3rd Qu.: 94.00
## Max. :760.00 Max. :1824.0 Max. :2997.0 Max. :1620.00
##
## HHVEH4 HHVEH5 HHVEH6 HHVEH7
## Min. : 0.00 Min. : 0.00 Min. : 0.000 Min. : 0.0000
## 1st Qu.: 10.00 1st Qu.: 1.00 1st Qu.: 0.000 1st Qu.: 0.0000
## Median : 17.00 Median : 3.00 Median : 1.000 Median : 0.0000
## Mean : 21.09 Mean : 4.95 Mean : 2.787 Mean : 0.7652
## 3rd Qu.: 27.00 3rd Qu.: 6.00 3rd Qu.: 3.000 3rd Qu.: 1.0000
## Max. :320.00 Max. :149.00 Max. :72.000 Max. :24.0000
##
## HHVEH8 notrav LPAgrp
## Min. : 0.0000 Min. : 0.0 Min. :0.000
## 1st Qu.: 0.0000 1st Qu.: 226.0 1st Qu.:3.000
## Median : 0.0000 Median : 321.0 Median :3.000
## Mean : 0.1916 Mean : 359.6 Mean :3.178
## 3rd Qu.: 0.0000 3rd Qu.: 447.0 3rd Qu.:4.000
## Max. :38.0000 Max. :3763.0 Max. :4.000
## Component 3: this is the bobox (bounding box of coordinates that is drawn around the boundaries of CA) Component 4: Is the proj4string that contains the projections.
The @ is used to access a specific slot of the spatial data frame
Operations on a column of data is allowed as usual
summary (CA.poly@data$VMT)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 21738 31289 37816 45815 1148168CA.poly@data$VMTpr = CA.poly@data$VMT/CA.poly@data$n_pr
summary (CA.poly@data$VMTpr)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00 18.49 22.46 24.06 28.37 90.00 63First define dummy variables that are based on the type of land use in each block group
CA.poly@data$center = CA.poly@data$LPAgrp == 4
CA.poly@data$suburb = CA.poly@data$LPAgrp == 3
CA.poly@data$exurb = CA.poly@data$LPAgrp == 2
CA.poly@data$rural = CA.poly@data$LPAgrp == 1
CA.poly@data$none = CA.poly@data$LPAgrp == 0Then estimate a regression model.
VMTprOLS<-lm(VMTpr~suburb+exurb+rural+HHVEH0+HHVEH1+HHVEH2+HHVEH3+HHVEH4+HHVEH5+HHVEH6+HHAGE7, data=CA.poly@data)
VMTprOLS##
## Call:
## lm(formula = VMTpr ~ suburb + exurb + rural + HHVEH0 + HHVEH1 +
## HHVEH2 + HHVEH3 + HHVEH4 + HHVEH5 + HHVEH6 + HHAGE7, data = CA.poly@data)
##
## Coefficients:
## (Intercept) suburbTRUE exurbTRUE ruralTRUE HHVEH0
## 19.474709 5.599604 5.391621 11.671735 -0.019978
## HHVEH1 HHVEH2 HHVEH3 HHVEH4 HHVEH5
## -0.009254 0.010847 0.010957 0.046618 -0.147416
## HHVEH6 HHAGE7
## -0.064597 0.002544summary(VMTprOLS)##
## Call:
## lm(formula = VMTpr ~ suburb + exurb + rural + HHVEH0 + HHVEH1 +
## HHVEH2 + HHVEH3 + HHVEH4 + HHVEH5 + HHVEH6 + HHAGE7, data = CA.poly@data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.705 -3.446 -0.754 2.503 63.486
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.4747095 0.1044064 186.528 < 2e-16 ***
## suburbTRUE 5.5996042 0.1005532 55.688 < 2e-16 ***
## exurbTRUE 5.3916213 0.1555270 34.667 < 2e-16 ***
## ruralTRUE 11.6717350 0.2011197 58.034 < 2e-16 ***
## HHVEH0 -0.0199781 0.0017925 -11.145 < 2e-16 ***
## HHVEH1 -0.0092536 0.0007676 -12.055 < 2e-16 ***
## HHVEH2 0.0108465 0.0008535 12.708 < 2e-16 ***
## HHVEH3 0.0109570 0.0019188 5.710 1.14e-08 ***
## HHVEH4 0.0466177 0.0040481 11.516 < 2e-16 ***
## HHVEH5 -0.1474161 0.0087433 -16.861 < 2e-16 ***
## HHVEH6 -0.0645967 0.0099456 -6.495 8.47e-11 ***
## HHAGE7 0.0025441 0.0011925 2.133 0.0329 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.687 on 23123 degrees of freedom
## (63 observations deleted due to missingness)
## Multiple R-squared: 0.4257, Adjusted R-squared: 0.4254
## F-statistic: 1558 on 11 and 23123 DF, p-value: < 2.2e-16Test the residuals from this regression model
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(sandwich) # I need this for the robust vcov matrixcoeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "const")) # this is the same as in least squares with no White adjsustment##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.47470947 0.10440640 186.5279 < 2.2e-16 ***
## suburbTRUE 5.59960420 0.10055324 55.6880 < 2.2e-16 ***
## exurbTRUE 5.39162132 0.15552699 34.6668 < 2.2e-16 ***
## ruralTRUE 11.67173496 0.20111973 58.0338 < 2.2e-16 ***
## HHVEH0 -0.01997811 0.00179251 -11.1453 < 2.2e-16 ***
## HHVEH1 -0.00925356 0.00076761 -12.0550 < 2.2e-16 ***
## HHVEH2 0.01084651 0.00085354 12.7077 < 2.2e-16 ***
## HHVEH3 0.01095696 0.00191877 5.7104 1.141e-08 ***
## HHVEH4 0.04661772 0.00404807 11.5160 < 2.2e-16 ***
## HHVEH5 -0.14741614 0.00874325 -16.8606 < 2.2e-16 ***
## HHVEH6 -0.06459674 0.00994559 -6.4950 8.472e-11 ***
## HHAGE7 0.00254406 0.00119249 2.1334 0.0329 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC0")) # this is the traditional White's adjustment to the var-Cov of the coefficient estimates##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.47470947 0.10555576 184.4969 < 2.2e-16 ***
## suburbTRUE 5.59960420 0.10105296 55.4126 < 2.2e-16 ***
## exurbTRUE 5.39162132 0.15242269 35.3728 < 2.2e-16 ***
## ruralTRUE 11.67173496 0.20380432 57.2693 < 2.2e-16 ***
## HHVEH0 -0.01997811 0.00209086 -9.5550 < 2.2e-16 ***
## HHVEH1 -0.00925356 0.00084993 -10.8874 < 2.2e-16 ***
## HHVEH2 0.01084651 0.00106403 10.1938 < 2.2e-16 ***
## HHVEH3 0.01095696 0.00211560 5.1791 2.248e-07 ***
## HHVEH4 0.04661772 0.00558563 8.3460 < 2.2e-16 ***
## HHVEH5 -0.14741614 0.01044866 -14.1086 < 2.2e-16 ***
## HHVEH6 -0.06459674 0.00796668 -8.1084 5.383e-16 ***
## HHAGE7 0.00254406 0.00118435 2.1481 0.03172 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# the following are other versions of computing the variance matrix of coefficient estimates
coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC1")) # this is improved White's adjustment ##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.47470947 0.10558315 184.4490 < 2.2e-16 ***
## suburbTRUE 5.59960420 0.10107918 55.3982 < 2.2e-16 ***
## exurbTRUE 5.39162132 0.15246224 35.3637 < 2.2e-16 ***
## ruralTRUE 11.67173496 0.20385720 57.2545 < 2.2e-16 ***
## HHVEH0 -0.01997811 0.00209140 -9.5525 < 2.2e-16 ***
## HHVEH1 -0.00925356 0.00085015 -10.8846 < 2.2e-16 ***
## HHVEH2 0.01084651 0.00106430 10.1912 < 2.2e-16 ***
## HHVEH3 0.01095696 0.00211615 5.1778 2.264e-07 ***
## HHVEH4 0.04661772 0.00558708 8.3438 < 2.2e-16 ***
## HHVEH5 -0.14741614 0.01045138 -14.1050 < 2.2e-16 ***
## HHVEH6 -0.06459674 0.00796875 -8.1063 5.477e-16 ***
## HHAGE7 0.00254406 0.00118465 2.1475 0.03176 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC2")) # this is another improved White's adjustment ##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.47470947 0.10582176 184.0331 < 2.2e-16 ***
## suburbTRUE 5.59960420 0.10120667 55.3284 < 2.2e-16 ***
## exurbTRUE 5.39162132 0.15263767 35.3230 < 2.2e-16 ***
## ruralTRUE 11.67173496 0.20403023 57.2059 < 2.2e-16 ***
## HHVEH0 -0.01997811 0.00209670 -9.5284 < 2.2e-16 ***
## HHVEH1 -0.00925356 0.00085103 -10.8734 < 2.2e-16 ***
## HHVEH2 0.01084651 0.00106576 10.1773 < 2.2e-16 ***
## HHVEH3 0.01095696 0.00212029 5.1677 2.390e-07 ***
## HHVEH4 0.04661772 0.00559303 8.3350 < 2.2e-16 ***
## HHVEH5 -0.14741614 0.01049748 -14.0430 < 2.2e-16 ***
## HHVEH6 -0.06459674 0.00798575 -8.0890 6.309e-16 ***
## HHAGE7 0.00254406 0.00118731 2.1427 0.03215 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC3")) # this is the third improvement##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.47470947 0.10609084 183.5664 < 2.2e-16 ***
## suburbTRUE 5.59960420 0.10136284 55.2432 < 2.2e-16 ***
## exurbTRUE 5.39162132 0.15285505 35.2728 < 2.2e-16 ***
## ruralTRUE 11.67173496 0.20425785 57.1422 < 2.2e-16 ***
## HHVEH0 -0.01997811 0.00210260 -9.5016 < 2.2e-16 ***
## HHVEH1 -0.00925356 0.00085214 -10.8592 < 2.2e-16 ***
## HHVEH2 0.01084651 0.00106751 10.1606 < 2.2e-16 ***
## HHVEH3 0.01095696 0.00212513 5.1559 2.545e-07 ***
## HHVEH4 0.04661772 0.00560063 8.3237 < 2.2e-16 ***
## HHVEH5 -0.14741614 0.01054746 -13.9765 < 2.2e-16 ***
## HHVEH6 -0.06459674 0.00800508 -8.0695 7.401e-16 ***
## HHAGE7 0.00254406 0.00119031 2.1373 0.03258 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(VMTprOLS, vcov = vcovHC(VMTprOLS, type = "HC4")) # this is the fourth improvement ##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.47470947 0.10660493 182.6811 < 2.2e-16 ***
## suburbTRUE 5.59960420 0.10164097 55.0920 < 2.2e-16 ***
## exurbTRUE 5.39162132 0.15320295 35.1927 < 2.2e-16 ***
## ruralTRUE 11.67173496 0.20456234 57.0571 < 2.2e-16 ***
## HHVEH0 -0.01997811 0.00211390 -9.4508 < 2.2e-16 ***
## HHVEH1 -0.00925356 0.00085395 -10.8362 < 2.2e-16 ***
## HHVEH2 0.01084651 0.00107056 10.1317 < 2.2e-16 ***
## HHVEH3 0.01095696 0.00213433 5.1337 2.864e-07 ***
## HHVEH4 0.04661772 0.00561357 8.3045 < 2.2e-16 ***
## HHVEH5 -0.14741614 0.01064813 -13.8443 < 2.2e-16 ***
## HHVEH6 -0.06459674 0.00804135 -8.0331 9.957e-16 ***
## HHAGE7 0.00254406 0.00119597 2.1272 0.03341 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1bptest(VMTprOLS, studentize=FALSE ) # the original BP test Page 62 class notes##
## Breusch-Pagan test
##
## data: VMTprOLS
## BP = 4989.3, df = 11, p-value < 2.2e-16bptest(VMTprOLS, studentize=TRUE) # the studentized version that is more robust to non-normal residuals##
## studentized Breusch-Pagan test
##
## data: VMTprOLS
## BP = 1092.9, df = 11, p-value < 2.2e-16If I run a typical serial correlation test, I may find that there is no correlation between pairs of rows We will check this later
dwtest(VMTprOLS)##
## Durbin-Watson test
##
## data: VMTprOLS
## DW = 1.9262, p-value = 9.589e-09
## alternative hypothesis: true autocorrelation is greater than 0Part 2 Extract Riverside County Data to Perform Analysis of Census Block Groups
We can select a portion of the data. For example, selecting two counties TWOCOUNTY <- CA.poly[CA.poly@data$countyname== c(“Riverside” , “Impreial”), ]
For the class example I just want to work with one county
YCOUNTY <- CA.poly[CA.poly@data$countyname== c("Riverside"), ]Then I want to know what is in the data slot of this county If you just run YCOUNTY@data you will get all the records - Try just for fun
You can give this instruction to look at the content of the data (not run here to keep output small)
summary(YCOUNTY@data)
Let’s estimate the same regression model we estimated for the entire State of California
VMTprOLSRiver<-lm(VMTpr~suburb+exurb+rural+HHVEH0+HHVEH1+HHVEH2+HHVEH3+HHVEH4+HHVEH5+HHVEH6+HHAGE7, data=YCOUNTY@data)
summary(VMTprOLSRiver)##
## Call:
## lm(formula = VMTpr ~ suburb + exurb + rural + HHVEH0 + HHVEH1 +
## HHVEH2 + HHVEH3 + HHVEH4 + HHVEH5 + HHVEH6 + HHAGE7, data = YCOUNTY@data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.319 -3.779 -0.688 3.021 49.862
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.687987 0.734600 29.524 < 2e-16 ***
## suburbTRUE 5.004393 0.688883 7.265 7.44e-13 ***
## exurbTRUE 2.831428 0.862265 3.284 0.00106 **
## ruralTRUE 8.050345 1.153020 6.982 5.25e-12 ***
## HHVEH0 -0.039912 0.009968 -4.004 6.68e-05 ***
## HHVEH1 -0.018997 0.003620 -5.248 1.87e-07 ***
## HHVEH2 0.001059 0.003655 0.290 0.77216
## HHVEH3 0.020531 0.007350 2.793 0.00532 **
## HHVEH4 0.112590 0.015355 7.332 4.61e-13 ***
## HHVEH5 -0.313633 0.038005 -8.253 4.77e-16 ***
## HHVEH6 0.107290 0.056601 1.896 0.05830 .
## HHAGE7 0.017379 0.003758 4.625 4.23e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.685 on 1017 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.4126, Adjusted R-squared: 0.4062
## F-statistic: 64.93 on 11 and 1017 DF, p-value: < 2.2e-16This model is different than what we got using the entire State (see the coefficient for exurbTRUE and compare it to suburbTRUE)
You can also do the BP and DW tests and see what you get but we can do some other more intersting things.
Building Neighborhoods Using Contiguity Rules. The example below uses queen moves (from Chess) to build links among block groups It also creates centroids as centers of “gravity”
In R this function is poly2nb
This builds a neighbours list based on regions with contiguous boundaries. queen=T allows even for a single polygon neighbor to meet the contiguity condition.
list.queenY<-poly2nb(YCOUNTY, queen=T)
coordsY<-coordinates(YCOUNTY)
plot(YCOUNTY)
plot(list.queenY, coordsY, add=T)
summary(list.queenY)## Neighbour list object:
## Number of regions: 1030
## Number of nonzero links: 6508
## Percentage nonzero weights: 0.6134414
## Average number of links: 6.318447
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 22
## 1 12 56 102 215 217 186 113 62 28 16 9 3 3 3 2 2
## 1 least connected region:
## 9993 with 1 link
## 2 most connected regions:
## 20444 20592 with 22 linksUsing neighborhoods based on the k-nearest neighbor rule
coords<-coordinates(YCOUNTY)
IDs<-row.names(as(YCOUNTY, "data.frame"))
plot(YCOUNTY)
sids_kn10<-knn2nb(knearneigh(coords, k=10), row.names=IDs)
plot(sids_kn10, coordsY, add=T)
summary (sids_kn10)## Neighbour list object:
## Number of regions: 1030
## Number of nonzero links: 10300
## Percentage nonzero weights: 0.9708738
## Average number of links: 10
## Non-symmetric neighbours list
## Link number distribution:
##
## 10
## 1030
## 1030 least connected regions:
## 9 101 107 134 174 318 364 420 454 493 501 520 525 533 546 563 575 628 642 670 695 699 724 742 821 822 886 1130 1198 1233 1239 1279 1290 1295 1312 1345 1400 1405 1411 1441 1454 1506 1527 1553 1556 1609 1650 1656 1691 1697 1703 1714 1744 1746 1790 1813 1844 1884 1900 1905 1981 1995 2033 2079 2104 2249 2253 2270 2289 2298 2326 2352 3964 3965 3966 3967 3968 3969 3970 3971 3972 3973 3974 3975 3976 3977 3978 3979 3980 3981 3982 3983 3984 3985 3986 3987 3988 3989 4001 4002 4003 4004 4005 4111 4112 4113 4114 4115 4116 4117 4118 4119 4120 4121 4122 4123 4124 4125 4126 4127 4128 4129 4130 4131 4132 4133 4134 4135 4136 4137 4138 4139 4140 4141 4142 4143 4144 4145 4146 4147 4148 4149 4150 4151 4702 4703 4704 4705 4706 4707 4708 4709 4710 4711 4712 4713 4714 4715 4716 4717 4718 4719 4720 4721 4722 4723 4724 4725 4726 4727 4728 4729 4730 4731 4732 4733 4734 4735 4736 4737 4738 4739 4740 4741 4742 4743 4744 4745 4746 4747 4748 4749 4750 4751 4752 4753 4754 4755 4756 4757 4758 4759 4760 4761 4762 4763 4764 4765 4766 4767 4768 4769 4770 4771 4772 4773 4774 4775 4776 4777 4778 4779 4780 5356 5412 5419 5443 5509 5516 5563 5609 5735 5739 5872 5911 5947 5992 6081 6082 6103 6116 6128 6142 6153 6194 6235 6253 6264 6268 6270 6278 6296 6344 6369 6379 6382 6428 6437 6457 6466 6467 6490 6492 6509 6513 6514 6546 6549 6587 6599 6603 6655 6702 6736 6791 6851 6864 6886 7024 7036 7176 7222 7260 7277 7295 7318 7400 7404 7510 7522 7801 8173 8174 8175 8176 8177 8178 8179 8184 8185 8186 8187 8188 8189 8190 8197 8198 8199 8200 8201 8202 8211 8212 8213 8214 8436 8437 8438 8439 8440 8441 8442 8443 8449 8450 8451 8452 8453 8454 8455 8456 8457 8458 8459 8462 8463 8466 8467 8468 8469 8470 8657 8660 8730 8747 8775 8800 8802 8896 8922 8930 8947 9032 9033 9035 9050 9099 9194 9199 9221 9224 9236 9257 9593 9607 9608 9609 9610 9623 9672 9673 9674 9679 9680 9681 9691 9692 9697 9698 9699 9700 9701 9702 9703 9704 9705 9706 9818 9866 9899 9932 9934 9993 10006 10026 10029 10033 10036 10049 10069 10100 10170 10180 10193 10198 10289 10345 10394 10490 10491 10492 10493 10495 10525 10526 10531 10532 10533 10534 10536 10537 10538 10653 10661 10663 10672 10690 10704 10714 10754 10775 10794 10838 10863 10898 10900 10906 10907 10908 10910 10924 10928 11013 11044 11046 11071 11098 11112 11146 11158 11170 11172 11217 11232 11236 11242 11290 11292 11297 11303 11329 11332 11336 11344 11365 11371 11379 11388 11393 11429 11445 11520 11522 11527 11531 11552 11553 11554 11562 11629 11638 11701 11717 11719 11722 11769 11792 11803 11804 11825 11839 11894 11917 11922 11966 11985 11992 12030 12045 12069 12073 12122 12133 12145 12149 12251 12267 12288 12302 12306 12312 12344 12365 12406 12534 12542 12572 12608 12642 12651 12663 12681 12690 12701 12706 12708 12766 12775 12780 12915 12916 12928 12947 12953 13004 13034 13037 13056 13091 13096 13120 13128 13140 13148 13174 13179 13277 13291 13315 13364 13371 13390 13394 13426 13428 13448 13511 13548 13565 13730 13735 15531 15532 15533 15534 15535 15536 15537 15538 15539 15547 15548 15549 15550 15551 15565 15566 15567 15568 15569 15570 15571 15572 15573 15574 15575 15576 15577 15589 15590 15591 15592 15593 15594 15595 15596 15617 15618 15635 15636 15637 15638 15639 15640 15688 15689 15690 15691 15692 15707 15708 15709 15710 15711 15712 15713 15714 15715 15761 15762 15763 15764 15765 15766 15767 15768 15769 15770 15771 15772 15773 15774 15775 15796 15797 15798 15799 15800 15801 15802 15803 15804 15805 15843 15844 15845 15846 15847 15848 15849 15850 15860 15861 15862 15863 15864 15886 15887 16511 16512 16534 16535 16536 16537 16538 16539 16540 16541 16542 16543 16544 16545 16546 16547 16548 16549 16550 16551 16552 16553 16554 16555 16556 16557 16558 16559 16560 16561 16562 16563 16564 16565 16566 16567 16568 16569 16570 16571 16572 16573 16574 16575 16576 16577 16578 16579 16580 16581 16582 16583 16584 16585 16586 16587 16588 16589 16590 16591 16592 16593 16594 16595 16596 16597 16598 16599 16600 16601 16602 16603 16604 16605 16606 16607 16608 16609 16610 16611 16612 16613 16614 16615 16616 16617 16618 16619 16620 16621 16622 16623 16624 16625 16626 16627 16628 17383 17411 17424 17429 17518 17521 17525 17535 17572 17587 17600 17621 17625 17633 17704 17705 17708 17720 17721 17728 17731 17732 17787 17790 17812 17817 17830 17848 17913 17916 17932 17943 18005 18007 18035 18045 18059 18065 18097 18114 18151 18170 18249 18260 18296 18307 18339 18371 18386 18413 18459 18484 18490 18502 18523 18532 18543 18556 18558 18618 18651 18674 18677 18697 18719 18776 18793 18858 18876 18909 18931 18944 18970 19044 19055 19056 19085 19099 19134 19142 19147 19172 19203 19208 19226 19247 19339 19353 19356 19372 19374 19392 19409 19411 19439 19447 19481 19520 19536 19581 19672 19698 19783 19786 19787 19791 19837 19958 20132 20357 20398 20415 20416 20417 20418 20419 20420 20421 20444 20516 20517 20518 20519 20520 20521 20544 20545 20546 20556 20566 20567 20568 20588 20592 20593 20594 20595 20656 20679 20680 20681 20686 20687 20688 20689 20690 20691 20726 20727 20728 20729 20730 20731 20732 20733 20734 20735 20739 20740 20741 20742 20743 20774 20775 20776 20777 20778 20779 20782 20794 20795 20796 20797 21011 21078 21079 21093 21176 21191 21277 21282 21294 21321 21326 21332 21359 21444 21447 21449 21463 21477 21509 21525 21530 21588 21607 21695 21727 21922 21926 21927 21928 21930 21931 21932 21935 21943 21944 21955 21970 21971 21972 21983 21984 22008 22022 22023 22024 22025 22026 22027 22028 22037 22038 22039 22040 22041 22042 22064 22083 22094 22177 22242 22300 22311 22323 22360 22428 22443 22450 22456 22462 22496 22505 22532 22564 22682 22683 22690 22743 22778 22815 22816 22817 22818 22825 22871 22872 22888 22889 22913 22923 22924 22925 22953 22988 23001 23011 23019 23107 23113 23116 23120 23164 23165 with 10 links
## 1030 most connected regions:
## 9 101 107 134 174 318 364 420 454 493 501 520 525 533 546 563 575 628 642 670 695 699 724 742 821 822 886 1130 1198 1233 1239 1279 1290 1295 1312 1345 1400 1405 1411 1441 1454 1506 1527 1553 1556 1609 1650 1656 1691 1697 1703 1714 1744 1746 1790 1813 1844 1884 1900 1905 1981 1995 2033 2079 2104 2249 2253 2270 2289 2298 2326 2352 3964 3965 3966 3967 3968 3969 3970 3971 3972 3973 3974 3975 3976 3977 3978 3979 3980 3981 3982 3983 3984 3985 3986 3987 3988 3989 4001 4002 4003 4004 4005 4111 4112 4113 4114 4115 4116 4117 4118 4119 4120 4121 4122 4123 4124 4125 4126 4127 4128 4129 4130 4131 4132 4133 4134 4135 4136 4137 4138 4139 4140 4141 4142 4143 4144 4145 4146 4147 4148 4149 4150 4151 4702 4703 4704 4705 4706 4707 4708 4709 4710 4711 4712 4713 4714 4715 4716 4717 4718 4719 4720 4721 4722 4723 4724 4725 4726 4727 4728 4729 4730 4731 4732 4733 4734 4735 4736 4737 4738 4739 4740 4741 4742 4743 4744 4745 4746 4747 4748 4749 4750 4751 4752 4753 4754 4755 4756 4757 4758 4759 4760 4761 4762 4763 4764 4765 4766 4767 4768 4769 4770 4771 4772 4773 4774 4775 4776 4777 4778 4779 4780 5356 5412 5419 5443 5509 5516 5563 5609 5735 5739 5872 5911 5947 5992 6081 6082 6103 6116 6128 6142 6153 6194 6235 6253 6264 6268 6270 6278 6296 6344 6369 6379 6382 6428 6437 6457 6466 6467 6490 6492 6509 6513 6514 6546 6549 6587 6599 6603 6655 6702 6736 6791 6851 6864 6886 7024 7036 7176 7222 7260 7277 7295 7318 7400 7404 7510 7522 7801 8173 8174 8175 8176 8177 8178 8179 8184 8185 8186 8187 8188 8189 8190 8197 8198 8199 8200 8201 8202 8211 8212 8213 8214 8436 8437 8438 8439 8440 8441 8442 8443 8449 8450 8451 8452 8453 8454 8455 8456 8457 8458 8459 8462 8463 8466 8467 8468 8469 8470 8657 8660 8730 8747 8775 8800 8802 8896 8922 8930 8947 9032 9033 9035 9050 9099 9194 9199 9221 9224 9236 9257 9593 9607 9608 9609 9610 9623 9672 9673 9674 9679 9680 9681 9691 9692 9697 9698 9699 9700 9701 9702 9703 9704 9705 9706 9818 9866 9899 9932 9934 9993 10006 10026 10029 10033 10036 10049 10069 10100 10170 10180 10193 10198 10289 10345 10394 10490 10491 10492 10493 10495 10525 10526 10531 10532 10533 10534 10536 10537 10538 10653 10661 10663 10672 10690 10704 10714 10754 10775 10794 10838 10863 10898 10900 10906 10907 10908 10910 10924 10928 11013 11044 11046 11071 11098 11112 11146 11158 11170 11172 11217 11232 11236 11242 11290 11292 11297 11303 11329 11332 11336 11344 11365 11371 11379 11388 11393 11429 11445 11520 11522 11527 11531 11552 11553 11554 11562 11629 11638 11701 11717 11719 11722 11769 11792 11803 11804 11825 11839 11894 11917 11922 11966 11985 11992 12030 12045 12069 12073 12122 12133 12145 12149 12251 12267 12288 12302 12306 12312 12344 12365 12406 12534 12542 12572 12608 12642 12651 12663 12681 12690 12701 12706 12708 12766 12775 12780 12915 12916 12928 12947 12953 13004 13034 13037 13056 13091 13096 13120 13128 13140 13148 13174 13179 13277 13291 13315 13364 13371 13390 13394 13426 13428 13448 13511 13548 13565 13730 13735 15531 15532 15533 15534 15535 15536 15537 15538 15539 15547 15548 15549 15550 15551 15565 15566 15567 15568 15569 15570 15571 15572 15573 15574 15575 15576 15577 15589 15590 15591 15592 15593 15594 15595 15596 15617 15618 15635 15636 15637 15638 15639 15640 15688 15689 15690 15691 15692 15707 15708 15709 15710 15711 15712 15713 15714 15715 15761 15762 15763 15764 15765 15766 15767 15768 15769 15770 15771 15772 15773 15774 15775 15796 15797 15798 15799 15800 15801 15802 15803 15804 15805 15843 15844 15845 15846 15847 15848 15849 15850 15860 15861 15862 15863 15864 15886 15887 16511 16512 16534 16535 16536 16537 16538 16539 16540 16541 16542 16543 16544 16545 16546 16547 16548 16549 16550 16551 16552 16553 16554 16555 16556 16557 16558 16559 16560 16561 16562 16563 16564 16565 16566 16567 16568 16569 16570 16571 16572 16573 16574 16575 16576 16577 16578 16579 16580 16581 16582 16583 16584 16585 16586 16587 16588 16589 16590 16591 16592 16593 16594 16595 16596 16597 16598 16599 16600 16601 16602 16603 16604 16605 16606 16607 16608 16609 16610 16611 16612 16613 16614 16615 16616 16617 16618 16619 16620 16621 16622 16623 16624 16625 16626 16627 16628 17383 17411 17424 17429 17518 17521 17525 17535 17572 17587 17600 17621 17625 17633 17704 17705 17708 17720 17721 17728 17731 17732 17787 17790 17812 17817 17830 17848 17913 17916 17932 17943 18005 18007 18035 18045 18059 18065 18097 18114 18151 18170 18249 18260 18296 18307 18339 18371 18386 18413 18459 18484 18490 18502 18523 18532 18543 18556 18558 18618 18651 18674 18677 18697 18719 18776 18793 18858 18876 18909 18931 18944 18970 19044 19055 19056 19085 19099 19134 19142 19147 19172 19203 19208 19226 19247 19339 19353 19356 19372 19374 19392 19409 19411 19439 19447 19481 19520 19536 19581 19672 19698 19783 19786 19787 19791 19837 19958 20132 20357 20398 20415 20416 20417 20418 20419 20420 20421 20444 20516 20517 20518 20519 20520 20521 20544 20545 20546 20556 20566 20567 20568 20588 20592 20593 20594 20595 20656 20679 20680 20681 20686 20687 20688 20689 20690 20691 20726 20727 20728 20729 20730 20731 20732 20733 20734 20735 20739 20740 20741 20742 20743 20774 20775 20776 20777 20778 20779 20782 20794 20795 20796 20797 21011 21078 21079 21093 21176 21191 21277 21282 21294 21321 21326 21332 21359 21444 21447 21449 21463 21477 21509 21525 21530 21588 21607 21695 21727 21922 21926 21927 21928 21930 21931 21932 21935 21943 21944 21955 21970 21971 21972 21983 21984 22008 22022 22023 22024 22025 22026 22027 22028 22037 22038 22039 22040 22041 22042 22064 22083 22094 22177 22242 22300 22311 22323 22360 22428 22443 22450 22456 22462 22496 22505 22532 22564 22682 22683 22690 22743 22778 22815 22816 22817 22818 22825 22871 22872 22888 22889 22913 22923 22924 22925 22953 22988 23001 23011 23019 23107 23113 23116 23120 23164 23165 with 10 linksHouse keeping task - Make sure we have valid data in all polygons
One complication is that one of the block groups has NA in VMTpr I change this to zero to avoid missing data and loose a polygon
summary(YCOUNTY@data$VMTpr)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 11.84 20.66 25.02 25.96 30.48 76.48 1YCOUNTY@data$VMTpr[is.na(YCOUNTY@data$VMTpr)] <- 0
summary(YCOUNTY@data$VMTpr)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 20.63 25.00 25.93 30.48 76.48Weights creation. This is the matrix W from lecture notes. We create two types of weights Queen and k=10 nearest neighbors weights
Weights without row standardization look like this
$weights[[670]][1] 1 1 1 1 1 1 1
$weights[[671]][1] 1 1 1 1 1 1 1
$weights[[672]][1] 1 1 1 1 1
$weights[[673]][1] 1 1 1 1 1
$weights[[674]][1] 1 1 1
$weights[[675]][1] 1 1 1 1 1 1
$weights[[676]][1] 1 1 1 1 1 1 1 1
Weights with row standardization look like this (I get this out using head(weightsname)
$weights[[679]][1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571
$weights[[680]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
$weights[[681]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
$weights[[682]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
$weights[[683]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
$weights[[684]][1] 0.25 0.25 0.25 0.25
$weights[[685]][1] 0.2 0.2 0.2 0.2 0.2
$weights[[686]][1] 0.2 0.2 0.2 0.2 0.2
$weights[[687]][1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
Create weights using 0 and 1s for connectivity.
From spdep: nb2listw(neighbours, glist=NULL, style=“codong type of weights”, zero.policy=FALSE)
The function adds a weights list with values given by the coding scheme style chosen. B is the basic binary coding, W is row standardised (sums over all links to n), C is globally standardised (sums over all links to n), U is equal to C divided by the number of neighbours (sums over all links to unity), while S is the variance-stabilizing coding scheme sums over all links to n.
queen_w <- nb2listw(list.queenY, style="B")
summary(queen_w)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 1030
## Number of nonzero links: 6508
## Percentage nonzero weights: 0.6134414
## Average number of links: 6.318447
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 22
## 1 12 56 102 215 217 186 113 62 28 16 9 3 3 3 2 2
## 1 least connected region:
## 9993 with 1 link
## 2 most connected regions:
## 20444 20592 with 22 links
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 1030 1060900 6508 13016 184088Create weights using row standardized weights
queen_ws <- nb2listw(list.queenY)
summary(queen_ws)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 1030
## Number of nonzero links: 6508
## Percentage nonzero weights: 0.6134414
## Average number of links: 6.318447
## Link number distribution:
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 22
## 1 12 56 102 215 217 186 113 62 28 16 9 3 3 3 2 2
## 1 least connected region:
## 9993 with 1 link
## 2 most connected regions:
## 20444 20592 with 22 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 1030 1060900 1030 345.3927 4291.3The same weight creation but using k nearest neighbor
sids_kn10_w<- nb2listw(sids_kn10)
summary(sids_kn10_w)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 1030
## Number of nonzero links: 10300
## Percentage nonzero weights: 0.9708738
## Average number of links: 10
## Non-symmetric neighbours list
## Link number distribution:
##
## 10
## 1030
## 1030 least connected regions:
## 9 101 107 134 174 318 364 420 454 493 501 520 525 533 546 563 575 628 642 670 695 699 724 742 821 822 886 1130 1198 1233 1239 1279 1290 1295 1312 1345 1400 1405 1411 1441 1454 1506 1527 1553 1556 1609 1650 1656 1691 1697 1703 1714 1744 1746 1790 1813 1844 1884 1900 1905 1981 1995 2033 2079 2104 2249 2253 2270 2289 2298 2326 2352 3964 3965 3966 3967 3968 3969 3970 3971 3972 3973 3974 3975 3976 3977 3978 3979 3980 3981 3982 3983 3984 3985 3986 3987 3988 3989 4001 4002 4003 4004 4005 4111 4112 4113 4114 4115 4116 4117 4118 4119 4120 4121 4122 4123 4124 4125 4126 4127 4128 4129 4130 4131 4132 4133 4134 4135 4136 4137 4138 4139 4140 4141 4142 4143 4144 4145 4146 4147 4148 4149 4150 4151 4702 4703 4704 4705 4706 4707 4708 4709 4710 4711 4712 4713 4714 4715 4716 4717 4718 4719 4720 4721 4722 4723 4724 4725 4726 4727 4728 4729 4730 4731 4732 4733 4734 4735 4736 4737 4738 4739 4740 4741 4742 4743 4744 4745 4746 4747 4748 4749 4750 4751 4752 4753 4754 4755 4756 4757 4758 4759 4760 4761 4762 4763 4764 4765 4766 4767 4768 4769 4770 4771 4772 4773 4774 4775 4776 4777 4778 4779 4780 5356 5412 5419 5443 5509 5516 5563 5609 5735 5739 5872 5911 5947 5992 6081 6082 6103 6116 6128 6142 6153 6194 6235 6253 6264 6268 6270 6278 6296 6344 6369 6379 6382 6428 6437 6457 6466 6467 6490 6492 6509 6513 6514 6546 6549 6587 6599 6603 6655 6702 6736 6791 6851 6864 6886 7024 7036 7176 7222 7260 7277 7295 7318 7400 7404 7510 7522 7801 8173 8174 8175 8176 8177 8178 8179 8184 8185 8186 8187 8188 8189 8190 8197 8198 8199 8200 8201 8202 8211 8212 8213 8214 8436 8437 8438 8439 8440 8441 8442 8443 8449 8450 8451 8452 8453 8454 8455 8456 8457 8458 8459 8462 8463 8466 8467 8468 8469 8470 8657 8660 8730 8747 8775 8800 8802 8896 8922 8930 8947 9032 9033 9035 9050 9099 9194 9199 9221 9224 9236 9257 9593 9607 9608 9609 9610 9623 9672 9673 9674 9679 9680 9681 9691 9692 9697 9698 9699 9700 9701 9702 9703 9704 9705 9706 9818 9866 9899 9932 9934 9993 10006 10026 10029 10033 10036 10049 10069 10100 10170 10180 10193 10198 10289 10345 10394 10490 10491 10492 10493 10495 10525 10526 10531 10532 10533 10534 10536 10537 10538 10653 10661 10663 10672 10690 10704 10714 10754 10775 10794 10838 10863 10898 10900 10906 10907 10908 10910 10924 10928 11013 11044 11046 11071 11098 11112 11146 11158 11170 11172 11217 11232 11236 11242 11290 11292 11297 11303 11329 11332 11336 11344 11365 11371 11379 11388 11393 11429 11445 11520 11522 11527 11531 11552 11553 11554 11562 11629 11638 11701 11717 11719 11722 11769 11792 11803 11804 11825 11839 11894 11917 11922 11966 11985 11992 12030 12045 12069 12073 12122 12133 12145 12149 12251 12267 12288 12302 12306 12312 12344 12365 12406 12534 12542 12572 12608 12642 12651 12663 12681 12690 12701 12706 12708 12766 12775 12780 12915 12916 12928 12947 12953 13004 13034 13037 13056 13091 13096 13120 13128 13140 13148 13174 13179 13277 13291 13315 13364 13371 13390 13394 13426 13428 13448 13511 13548 13565 13730 13735 15531 15532 15533 15534 15535 15536 15537 15538 15539 15547 15548 15549 15550 15551 15565 15566 15567 15568 15569 15570 15571 15572 15573 15574 15575 15576 15577 15589 15590 15591 15592 15593 15594 15595 15596 15617 15618 15635 15636 15637 15638 15639 15640 15688 15689 15690 15691 15692 15707 15708 15709 15710 15711 15712 15713 15714 15715 15761 15762 15763 15764 15765 15766 15767 15768 15769 15770 15771 15772 15773 15774 15775 15796 15797 15798 15799 15800 15801 15802 15803 15804 15805 15843 15844 15845 15846 15847 15848 15849 15850 15860 15861 15862 15863 15864 15886 15887 16511 16512 16534 16535 16536 16537 16538 16539 16540 16541 16542 16543 16544 16545 16546 16547 16548 16549 16550 16551 16552 16553 16554 16555 16556 16557 16558 16559 16560 16561 16562 16563 16564 16565 16566 16567 16568 16569 16570 16571 16572 16573 16574 16575 16576 16577 16578 16579 16580 16581 16582 16583 16584 16585 16586 16587 16588 16589 16590 16591 16592 16593 16594 16595 16596 16597 16598 16599 16600 16601 16602 16603 16604 16605 16606 16607 16608 16609 16610 16611 16612 16613 16614 16615 16616 16617 16618 16619 16620 16621 16622 16623 16624 16625 16626 16627 16628 17383 17411 17424 17429 17518 17521 17525 17535 17572 17587 17600 17621 17625 17633 17704 17705 17708 17720 17721 17728 17731 17732 17787 17790 17812 17817 17830 17848 17913 17916 17932 17943 18005 18007 18035 18045 18059 18065 18097 18114 18151 18170 18249 18260 18296 18307 18339 18371 18386 18413 18459 18484 18490 18502 18523 18532 18543 18556 18558 18618 18651 18674 18677 18697 18719 18776 18793 18858 18876 18909 18931 18944 18970 19044 19055 19056 19085 19099 19134 19142 19147 19172 19203 19208 19226 19247 19339 19353 19356 19372 19374 19392 19409 19411 19439 19447 19481 19520 19536 19581 19672 19698 19783 19786 19787 19791 19837 19958 20132 20357 20398 20415 20416 20417 20418 20419 20420 20421 20444 20516 20517 20518 20519 20520 20521 20544 20545 20546 20556 20566 20567 20568 20588 20592 20593 20594 20595 20656 20679 20680 20681 20686 20687 20688 20689 20690 20691 20726 20727 20728 20729 20730 20731 20732 20733 20734 20735 20739 20740 20741 20742 20743 20774 20775 20776 20777 20778 20779 20782 20794 20795 20796 20797 21011 21078 21079 21093 21176 21191 21277 21282 21294 21321 21326 21332 21359 21444 21447 21449 21463 21477 21509 21525 21530 21588 21607 21695 21727 21922 21926 21927 21928 21930 21931 21932 21935 21943 21944 21955 21970 21971 21972 21983 21984 22008 22022 22023 22024 22025 22026 22027 22028 22037 22038 22039 22040 22041 22042 22064 22083 22094 22177 22242 22300 22311 22323 22360 22428 22443 22450 22456 22462 22496 22505 22532 22564 22682 22683 22690 22743 22778 22815 22816 22817 22818 22825 22871 22872 22888 22889 22913 22923 22924 22925 22953 22988 23001 23011 23019 23107 23113 23116 23120 23164 23165 with 10 links
## 1030 most connected regions:
## 9 101 107 134 174 318 364 420 454 493 501 520 525 533 546 563 575 628 642 670 695 699 724 742 821 822 886 1130 1198 1233 1239 1279 1290 1295 1312 1345 1400 1405 1411 1441 1454 1506 1527 1553 1556 1609 1650 1656 1691 1697 1703 1714 1744 1746 1790 1813 1844 1884 1900 1905 1981 1995 2033 2079 2104 2249 2253 2270 2289 2298 2326 2352 3964 3965 3966 3967 3968 3969 3970 3971 3972 3973 3974 3975 3976 3977 3978 3979 3980 3981 3982 3983 3984 3985 3986 3987 3988 3989 4001 4002 4003 4004 4005 4111 4112 4113 4114 4115 4116 4117 4118 4119 4120 4121 4122 4123 4124 4125 4126 4127 4128 4129 4130 4131 4132 4133 4134 4135 4136 4137 4138 4139 4140 4141 4142 4143 4144 4145 4146 4147 4148 4149 4150 4151 4702 4703 4704 4705 4706 4707 4708 4709 4710 4711 4712 4713 4714 4715 4716 4717 4718 4719 4720 4721 4722 4723 4724 4725 4726 4727 4728 4729 4730 4731 4732 4733 4734 4735 4736 4737 4738 4739 4740 4741 4742 4743 4744 4745 4746 4747 4748 4749 4750 4751 4752 4753 4754 4755 4756 4757 4758 4759 4760 4761 4762 4763 4764 4765 4766 4767 4768 4769 4770 4771 4772 4773 4774 4775 4776 4777 4778 4779 4780 5356 5412 5419 5443 5509 5516 5563 5609 5735 5739 5872 5911 5947 5992 6081 6082 6103 6116 6128 6142 6153 6194 6235 6253 6264 6268 6270 6278 6296 6344 6369 6379 6382 6428 6437 6457 6466 6467 6490 6492 6509 6513 6514 6546 6549 6587 6599 6603 6655 6702 6736 6791 6851 6864 6886 7024 7036 7176 7222 7260 7277 7295 7318 7400 7404 7510 7522 7801 8173 8174 8175 8176 8177 8178 8179 8184 8185 8186 8187 8188 8189 8190 8197 8198 8199 8200 8201 8202 8211 8212 8213 8214 8436 8437 8438 8439 8440 8441 8442 8443 8449 8450 8451 8452 8453 8454 8455 8456 8457 8458 8459 8462 8463 8466 8467 8468 8469 8470 8657 8660 8730 8747 8775 8800 8802 8896 8922 8930 8947 9032 9033 9035 9050 9099 9194 9199 9221 9224 9236 9257 9593 9607 9608 9609 9610 9623 9672 9673 9674 9679 9680 9681 9691 9692 9697 9698 9699 9700 9701 9702 9703 9704 9705 9706 9818 9866 9899 9932 9934 9993 10006 10026 10029 10033 10036 10049 10069 10100 10170 10180 10193 10198 10289 10345 10394 10490 10491 10492 10493 10495 10525 10526 10531 10532 10533 10534 10536 10537 10538 10653 10661 10663 10672 10690 10704 10714 10754 10775 10794 10838 10863 10898 10900 10906 10907 10908 10910 10924 10928 11013 11044 11046 11071 11098 11112 11146 11158 11170 11172 11217 11232 11236 11242 11290 11292 11297 11303 11329 11332 11336 11344 11365 11371 11379 11388 11393 11429 11445 11520 11522 11527 11531 11552 11553 11554 11562 11629 11638 11701 11717 11719 11722 11769 11792 11803 11804 11825 11839 11894 11917 11922 11966 11985 11992 12030 12045 12069 12073 12122 12133 12145 12149 12251 12267 12288 12302 12306 12312 12344 12365 12406 12534 12542 12572 12608 12642 12651 12663 12681 12690 12701 12706 12708 12766 12775 12780 12915 12916 12928 12947 12953 13004 13034 13037 13056 13091 13096 13120 13128 13140 13148 13174 13179 13277 13291 13315 13364 13371 13390 13394 13426 13428 13448 13511 13548 13565 13730 13735 15531 15532 15533 15534 15535 15536 15537 15538 15539 15547 15548 15549 15550 15551 15565 15566 15567 15568 15569 15570 15571 15572 15573 15574 15575 15576 15577 15589 15590 15591 15592 15593 15594 15595 15596 15617 15618 15635 15636 15637 15638 15639 15640 15688 15689 15690 15691 15692 15707 15708 15709 15710 15711 15712 15713 15714 15715 15761 15762 15763 15764 15765 15766 15767 15768 15769 15770 15771 15772 15773 15774 15775 15796 15797 15798 15799 15800 15801 15802 15803 15804 15805 15843 15844 15845 15846 15847 15848 15849 15850 15860 15861 15862 15863 15864 15886 15887 16511 16512 16534 16535 16536 16537 16538 16539 16540 16541 16542 16543 16544 16545 16546 16547 16548 16549 16550 16551 16552 16553 16554 16555 16556 16557 16558 16559 16560 16561 16562 16563 16564 16565 16566 16567 16568 16569 16570 16571 16572 16573 16574 16575 16576 16577 16578 16579 16580 16581 16582 16583 16584 16585 16586 16587 16588 16589 16590 16591 16592 16593 16594 16595 16596 16597 16598 16599 16600 16601 16602 16603 16604 16605 16606 16607 16608 16609 16610 16611 16612 16613 16614 16615 16616 16617 16618 16619 16620 16621 16622 16623 16624 16625 16626 16627 16628 17383 17411 17424 17429 17518 17521 17525 17535 17572 17587 17600 17621 17625 17633 17704 17705 17708 17720 17721 17728 17731 17732 17787 17790 17812 17817 17830 17848 17913 17916 17932 17943 18005 18007 18035 18045 18059 18065 18097 18114 18151 18170 18249 18260 18296 18307 18339 18371 18386 18413 18459 18484 18490 18502 18523 18532 18543 18556 18558 18618 18651 18674 18677 18697 18719 18776 18793 18858 18876 18909 18931 18944 18970 19044 19055 19056 19085 19099 19134 19142 19147 19172 19203 19208 19226 19247 19339 19353 19356 19372 19374 19392 19409 19411 19439 19447 19481 19520 19536 19581 19672 19698 19783 19786 19787 19791 19837 19958 20132 20357 20398 20415 20416 20417 20418 20419 20420 20421 20444 20516 20517 20518 20519 20520 20521 20544 20545 20546 20556 20566 20567 20568 20588 20592 20593 20594 20595 20656 20679 20680 20681 20686 20687 20688 20689 20690 20691 20726 20727 20728 20729 20730 20731 20732 20733 20734 20735 20739 20740 20741 20742 20743 20774 20775 20776 20777 20778 20779 20782 20794 20795 20796 20797 21011 21078 21079 21093 21176 21191 21277 21282 21294 21321 21326 21332 21359 21444 21447 21449 21463 21477 21509 21525 21530 21588 21607 21695 21727 21922 21926 21927 21928 21930 21931 21932 21935 21943 21944 21955 21970 21971 21972 21983 21984 22008 22022 22023 22024 22025 22026 22027 22028 22037 22038 22039 22040 22041 22042 22064 22083 22094 22177 22242 22300 22311 22323 22360 22428 22443 22450 22456 22462 22496 22505 22532 22564 22682 22683 22690 22743 22778 22815 22816 22817 22818 22825 22871 22872 22888 22889 22913 22923 22924 22925 22953 22988 23001 23011 23019 23107 23113 23116 23120 23164 23165 with 10 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 1030 1060900 1030 180.98 4232.62I you use style=“C” gives equal weights to all connections and produces the following in terms of weights
$weights[[671]][1] 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667
$weights[[672]][1] 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667
$weights[[673]][1] 0.1582667 0.1582667 0.1582667 0.1582667 0.1582667
$weights[[674]][1] 0.1582667 0.1582667 0.1582667
Autocorrelations at different lags using the defaults in sp.correlogram.
The question we ask is: Do we see spatial correlation in the vehicle miles per person produce by each Censusblock groups? How far is this correlation extending? What difference does it make to use weights using Queen continguity vs k-nearest neighbor continguity?
The following are the plots from slides 57 and 58 of class notes
mor10q <- sp.correlogram(list.queenY, var=YCOUNTY@data$VMTpr, order=10, method="I")
plot(mor10q, main = "Moran's I with Queen Contiguity and Row Standardization")
mor10k <- sp.correlogram(sids_kn10, var=YCOUNTY@data$VMTpr, order=10, method="I", zero.policy=TRUE) # need zero policy because some polygons are not connected
plot(mor10k, main = "Moran I with knn10 Contiguity and Row Standardization")
ger10q <- sp.correlogram(list.queenY, var=YCOUNTY@data$VMTpr, order=10, method="C")
plot(ger10q, main = "Geary's C with Queen Contiguity and Row Standardization")
ger10k <- sp.correlogram(sids_kn10, var=YCOUNTY@data$VMTpr, order=10, method="C", zero.policy=TRUE) # need zero policy because some polygons are not connected
plot(ger10k, main = "Geary's C with knn10 Contiguity and Row Standardization")
We can also create statstical tests for spatial correlation based on the z-scores
The first uses a weights matrix that is derived from Queen continguity using the default spatial weights row standardization and computing the variance with analytical randomization
moran.test(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY))##
## Moran I test under randomisation
##
## data: YCOUNTY@data$VMTpr
## weights: nb2listw(list.queenY)
##
## Moran I statistic standard deviate = 14.137, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2529411655 -0.0009718173 0.0003225841In the Table above
standarddeviate=((0.2529411655-(-0.0009718173))/sqrt(0.0003225841))
standarddeviate## [1] 14.1372This indicates that we have a strong spatial correlation
We repeat the same but this time we assume normality
moran.test(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY), randomisation=FALSE)##
## Moran I test under normality
##
## data: YCOUNTY@data$VMTpr
## weights: nb2listw(list.queenY)
##
## Moran I statistic standard deviate = 14.117, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2529411655 -0.0009718173 0.0003235222The standard deviate is not very different from the previous result.
This usually happens when we have many units (the polygons)
The next question is: are there observations that have very high spatial correlations?
spdep has a function called moran.plot. This produces an object with some useful quantities For example it runs a regression of x (the variable we analyze) on wx (the weighted values of the neighbors of x) It also plots each x and wx pairs and the mean of x and wx
Below we use Queen continuity and style=“C” This weights up observations with many neighbors because all connections to all polygons take the same value. Polygons with 6 connections will be influneced by 2 more neighbors than polygons with 4 connections (all weighted with the same value (0.15… we saw before))
msp <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="C"), quiet=TRUE)
title("Moran scatterplot Riverside")
The following code is from Bivand/Pebesma/Gomez-Rubio and extracts from object msp the influence of each pair x and wx in the regression of wx on x. This is one way to identify which polygons influences more the global Moran’s index
It also plots a map
infl <- apply(msp$is.inf, 1, any)
x <- YCOUNTY@data$VMTpr
lhx <- cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=TRUE)
wx <- lag(nb2listw(list.queenY, style="C"), YCOUNTY@data$VMTpr)
lhwx <- cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=TRUE)
lhlh <- interaction(lhx, lhwx, infl, drop=TRUE)
cols <- rep(1, length(lhlh))
cols[lhlh == "None"] <- 1
cols[lhlh == "H.L.TRUE"] <- 2
cols[lhlh == "L.H.TRUE"] <- 3
cols[lhlh == "H.H.TRUE"] <- 4
plot(YCOUNTY, col=brewer.pal(4, "Accent")[cols])
# RSB quietening greys
legend("topright", legend=c("None", "HL", "LH", "HH"), fill=brewer.pal(4, "Accent"), bty="n", cex=0.8, y.intersp=0.8)
title("Block groups with influence")
Not very nice and cannot even tell where the polygons (blockgroups). Leaflet will do the job.
First I define the colors and labels for the categories of influence and then build a map
colsF <- factor(cols, labels=c("None", "HighLow", "LowHigh", "HighHigh"))
LHpal <- colorFactor(topo.colors(4), colsF)
popup <- paste0("<strong> BLOCKGROUP </strong>", IDs)
leaflet(YCOUNTY) %>%
addPolygons(stroke = FALSE, smoothFactor = 0.2, fillOpacity = 0.7, popup=popup,
color = ~LHpal(colsF)) %>% addTiles() %>% addLegend("topright", pal = LHpal, values = ~colsF,
title = "Influence",
opacity = 1)The definiton weights has a substantial impact on the influence of neighbors and ultimately on Moran’s I From spdep: B is the basic binary coding, W is row standardised (sums over all links to n), C is globally standardised (sums over all links to n), U is equal to C divided by the number of neighbours (sums over all links to unity), while S is the variance-stabilizing coding scheme sums over all links to n.
mspdefault <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY), quiet=TRUE)
title("Moran scatterplot Riverside Default Style Default")
mspW <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="W"), quiet=TRUE)
title("Moran scatterplot Riverside Style W")
mspC <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="C"), quiet=TRUE)
title("Moran scatterplot Riverside Style C")
mspB <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="B"), quiet=TRUE)
title("Moran scatterplot Riverside Style B")
mspU <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="U"), quiet=TRUE)
title("Moran scatterplot Riverside Style U")
mspS <- moran.plot(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="S"), quiet=TRUE)
title("Moran scatterplot Riverside Style S")
The equation for local Moran’s I (one of the Anselin LISA indicators) - see Gauchospace paper is below. In this book Bivan et al show the denominator division by n not n-1. Numerically does not matter. Division by n is the defaul in spdep
\[ I_i = \frac{(x_i-\bar{x})}{{∑_{k=1}^{n}(x_k-\bar{x})^2}/(n-1)}{∑_{j=1}^{n}w_{ij}(x_j-\bar{x})} \]
We get as many of these indicators as the number of units (blockgroups in our example)
localM1 <- as.data.frame(localmoran(YCOUNTY@data$VMTpr, listw=nb2listw(list.queenY, style="C")))
summary(localM1)## Ii E.Ii Var.Ii
## Min. :-2.26283 Min. :-0.0033837 Min. :0.02493
## 1st Qu.:-0.02368 1st Qu.:-0.0010766 1st Qu.:0.12415
## Median : 0.09788 Median :-0.0009228 Median :0.14884
## Mean : 0.24063 Mean :-0.0009718 Mean :0.15658
## 3rd Qu.: 0.42794 3rd Qu.:-0.0007690 3rd Qu.:0.17348
## Max. : 4.69817 Max. :-0.0001538 Max. :0.53725
## Z.Ii Pr(z > 0)
## Min. :-5.43029 Min. :0.0000
## 1st Qu.:-0.05803 1st Qu.:0.1280
## Median : 0.24817 Median :0.4020
## Mean : 0.61887 Mean :0.3710
## 3rd Qu.: 1.13600 3rd Qu.:0.5231
## Max. :11.78184 Max. :1.0000Store the Local Moran I and its zscores in the database
YCOUNTY@data$localM1 <- localM1[,1]
summary(YCOUNTY@data$localM1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.26283 -0.02368 0.09788 0.24063 0.42794 4.69817YCOUNTY@data$zscoreM1 <- localM1[,4]
summary(YCOUNTY@data$zscoreM1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -5.43029 -0.05803 0.24817 0.61887 1.13600 11.78184Map the observed x
qpalreds <- colorNumeric(
palette = "Reds",
domain = CA.poly@data$VMTpr)
leaflet(YCOUNTY) %>%
addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2,
color = ~qpalreds(VMTpr)) %>% addTiles() %>% addLegend("topright", pal = qpalreds, values = ~VMTpr,
title = "Observed Miles per Person",
labFormat = labelFormat(prefix = "Observed Miles per person day "),
opacity = 0.8)qpalM <- colorNumeric(
palette = "Blues",
domain = CA.poly@data$localM1)
leaflet(YCOUNTY) %>%
addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2,
color = ~qpalM(localM1)) %>% addTiles() %>% addLegend("topright", pal = qpalM, values = ~localM1,
title = "Local Moran I",
labFormat = labelFormat(prefix = "Local Moran I "),
opacity = 1.0)qpalZ <- colorNumeric(
palette = "Greens",
domain = CA.poly@data$zscoreM1)
leaflet(YCOUNTY) %>%
addPolygons(stroke = FALSE, fillOpacity = .8, smoothFactor = 0.2,
color = ~qpalZ(zscoreM1)) %>% addTiles() %>% addLegend("topright", pal = qpalZ, values = ~zscoreM1,
title = "Local Moran I Z score",
labFormat = labelFormat(prefix = "Local Moran I Z score "),
opacity = 0.9)